---
title: "nonhhcc_nb_censored_loaded_weighted model"
output:
  html_document:
    df_print: paged
---

```{r setup, include=FALSE}

model.name <- 'nonhhcc_nb_censored_loaded_weighted'

library(flexdashboard)
library(tidyverse)
library(here)
library(ggthemes)
library(cowplot)
library(broom)
library(broom.mixed)
library(brms)
library(tidybayes)
library(tictoc)

library(gt)
library(gtsummary)
library(kableExtra)

theme_set(theme_cowplot())
```

Create an output directory for the model, if it doesn't already exist

```{r}
out.dir <- here('bics-paper-code', 'out')
dir.create(file.path(out.dir, model.name))

out.dir <- file.path(out.dir, model.name)
```

```{r}
# use 4 cores to fit 4 chains in parallel
options(mc.cores=4)
```

```{r read-data}
df_formod <- readRDS(file=file.path(out.dir, '..', 'waves0to2-for-models.rds'))

source(here('bics-paper-code', 'code', 'model-coef-plot-helpers.R'))
```


# Fit the model to all of the data with city dummies

```{r fit-model}
tic(glue::glue("Fitting model {model.name}"))
nb_brms <- brm(
                  formula = num_cc_nonhh | cens(is_topcoded_cc_nonhh) + weights(weight_pooled) ~ 
                    agecat_w0*gender + 
                    w_hhsize +
                    race_ethnicity*wave +
                    city*wave +
                    reference_weekday,
                  family = 'negbinomial',
                  # constant prior for Philadelphia in Wave 0 (we have no data here)
                  prior = set_prior('constant(0)', coef='cityPhiladelphia'),
                  chains=4,
                  inits="0",
                  data=df_formod)
                 
toc()
```

```{r}
summary(nb_brms)
```
```{r}
plot(conditional_effects(nb_brms, "agecat_w0:gender"))
```
```{r}
plot(conditional_effects(nb_brms, "wave"))
```
```{r}
plot(conditional_effects(nb_brms, "w_hhsize"))
```

```{r}
ce_weekday <- conditional_effects(nb_brms, "reference_weekday")
plot(ce_weekday)
```


```{r}
plot(conditional_effects(nb_brms, "city:wave"))
```

```{r}
#ce_wave_raceeth <- conditional_effects(nb_brms, "wave:race_ethnicity")
ce_wave_raceeth <- conditional_effects(nb_brms, "race_ethnicity:wave")

# don't want to plot the 'unknown' coefficient, which is based on very few obs, so has large CIs
# and is not interesting to us
ce_wave_raceeth_toplot <- ce_wave_raceeth
ce_wave_raceeth_toplot[[1]] <- ce_wave_raceeth_toplot[[1]] %>% 
  filter(effect2__ != "(Unknown)") %>%
  filter(effect1__ != "(Unknown)")

#plot(conditional_effects(nb_brms, "race_ethnicity:wave"))
plot(ce_wave_raceeth_toplot)
```

```{r}
ce_wave_raceeth <- conditional_effects(nb_brms, "wave:race_ethnicity")
#ce_wave_raceeth <- conditional_effects(nb_brms, "race_ethnicity:wave")

# don't want to plot the 'unknown' coefficient, which is based on very few obs, so has large CIs
# and is not interesting to us
ce_wave_raceeth_toplot <- ce_wave_raceeth
ce_wave_raceeth_toplot[[1]] <- ce_wave_raceeth_toplot[[1]] %>% 
  filter(effect2__ != "(Unknown)") %>%
  filter(effect1__ != "(Unknown)")

#plot(conditional_effects(nb_brms, "race_ethnicity:wave"))
plot(ce_wave_raceeth_toplot)
```

## Assess fit

Posterior predictive check

```{r}
pp_check(nb_brms, 
         type='hist')
```


## Inferences

```{r}
get_variables(nb_brms)
```

```{r plot-censoring-allcc, fig.height=10, fig.width=10}
m1_coef_plot <- plot_nb_res(nb_brms)$coef_plot

m1_coef_plot
```


```{r}
saveRDS(nb_brms, file.path(out.dir, paste0(model.name, '.rds')))
```

## Technical summary

### Priors

```{r}
prior_summary(nb_brms)
```

### Stan code generated by brms

```{r}
stancode(nb_brms)
```
